
    ##############################################################################################
    # R code to replicate
    # Figure 1 of Measuring Foreign Policy Positions of Members of the US Congress
    ##############################################################################################

    Bills <- read.csv("FP_All_H103.csv", header=T) # Foreign Policy Roll Call Data for the 103rd House
    Others <- read.csv("House103_Other_Measures.csv", header=T) # Other measures of ideology for the 103rd House

    library(MCMCpack)

    # prepare the roll call data for MCMCirt1d
    rownames(Bills) <- paste("A", Bills$ICPSR_ID, sep="")

    Sen  <- Bills[, -c(1:5)] # vote data
    N <- nrow(Sen)
    M <- ncol(Sen) 
    row.names(Sen) <- row.names(Bills)    
                                         
    MCMC <- MCMCirt1d(Sen, theta.constraints=list(A13035 = "-"),  # Charles Rangel (D-NY)
                      mcmc=10000, burnin=1000, thin=10, store.item=TRUE, store.ability=TRUE)   

    iter <- nrow(MCMC)
    
    ### Pick MCMC chains for ideal points and get the summary
    x.mc      <- MCMC[ ,1:N]      
    x.mean    <- round(apply(x.mc, 2, FUN=mean), digits=3)
    x.CI      <- round(apply(x.mc, 2, FUN=quantile, probs=c(0.025,0.975)), digits=3)
    IRT_Summary <- cbind(Bills[,1:5], x.mean, t(x.CI))
    colnames(IRT_Summary) <- c("ICPSR_ID", "Name", "ICPSR_state", "State_name", "Party", "IRT_mean", "IRT_2.5Q", "IRT_97.5Q")
    
    
    # Merge with other measures of ideologies
    All <- merge(Others, IRT_Summary, by="ICPSR_ID", all=F)

    ##### plot 
     
    postscript("Figure1.eps", height=8.5, width=8, horizontal=F)
    par(mfrow=c(2, 2))
    plot(All$IRT_mean, All$W_NOM1, type='p', col='black', xlim=c(-3,3), ylim=c(-1, 1), xlab="IRT Estimate", ylab="W-NOMINATE(dim 1)")
    segments(-3, -1, 3, 1, lty=1)
    plot(All$IRT_mean, -All$ADA_score, type='p', col='black', xlim=c(-3,3), ylim=c(min(-All$ADA_score), max(-All$ADA_score)), xlab="IRT Estimate", ylab="ADA Score (reversed)")
    segments(-3, -100, 3, 0, lty=1)
    plot(All$IRT_mean, All$ASC_score, type='p', col='black', xlim=c(-3,3), ylim=c(0, 100), xlab="IRT Estimate", ylab="ACS Score")
    segments(-3, 0, 3, 100, lty=1)
    plot(All$IRT_mean, All$NJ_FP, type='p', col='black', xlim=c(-3,3), ylim=c(0, 100), xlab="IRT Estimate", ylab="NJ Foreign Policy Score")
    segments(-3, 0, 3, 100, lty=1)
    dev.off()
    
